home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / umich / apps / math / lpsolves.zoo / solve.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-08-07  |  19.3 KB  |  876 lines
  1. #include "defines.h"
  2. #include "globals.h"
  3.  
  4.  
  5. void  condensecol(int rownr,
  6.           int numc,
  7.           double *pcol)
  8. {
  9.   int  i, elnr;
  10.   
  11.   if (Verbose)
  12.     printf("condensecol\n");
  13.   elnr = Endetacol[numc];
  14.  
  15.   if (elnr + Rows + 2 > Cur_eta_size) /* maximum local growth of Eta */
  16.     resize_eta();
  17.  
  18.   for (i = 0; i <= Rows; i++)
  19.     if (i != rownr && pcol[i] != 0)
  20.       {
  21.     Eta_rownr[elnr] = i;
  22.     Eta_value[elnr] = pcol[i];
  23.     elnr++;
  24.       }
  25.   Eta_rownr[elnr] = rownr;
  26.   Eta_value[elnr] = pcol[rownr];
  27.   elnr++;
  28.   Endetacol[numc + 1] = elnr;
  29. } /* condensecol */
  30.  
  31. void  addetacol(int *numc)
  32. {
  33.   int    i, j, k;
  34.   double theta;
  35.   
  36.   if (Verbose)
  37.     printf("addetacol\n");
  38.   j = Endetacol[(*numc)];
  39.   (*numc)++;
  40.   k = Endetacol[(*numc)];
  41.   theta = 1 / (double) Eta_value[k - 1];
  42.   Eta_value[k - 1] = theta;
  43.   for (i = j; i < k - 1; i++)
  44.     Eta_value[i] *= -theta;
  45. } /* addetacol */
  46.  
  47.  
  48. void  presolve(void)
  49. {
  50.   int    i, j, rownr;
  51.   double theta;
  52.   int    *num, *rownum;
  53.  
  54.   CALLOC(num, Rows + 2, int);
  55.   CALLOC(rownum, Rows + 2, int);
  56.   
  57.   if (Verbose)
  58.     printf("presolve\n");
  59.   Rh[0] = 0;
  60.   for (i = 0; i <= Sum; i++)
  61.     Lower[i] = TRUE;
  62.   for (i = 1; i <= Rows; i++)
  63.     Basis[i] = TRUE;
  64.   for (i = 1; i <= Columns; i++)
  65.     Basis[Rows + i] = FALSE;
  66.   for (i = 0; i <= Rows; i++)
  67.     Bas[i] = i;
  68.  
  69.   /* we want to maximise the objective function, so: */
  70.   Chsign[0] = TRUE;
  71.  
  72.   for (i = 1; i <= Rows; i++)
  73.     if ((Relat[i] == GE) && (Upbo[i] == INFINITE))
  74.       Chsign[i] = TRUE;
  75.     else
  76.       Chsign[i] = FALSE;
  77.  
  78.   /* invert values in rows with Chsign[i] == TRUE */
  79.   for (i = 1; i <= Rows; i++)
  80.     if (Chsign[i])
  81.       Rh[i] = -Rh[i];
  82.   for (i = 0; i < Nonnuls; i++)
  83.     if (Chsign[Mat[i].rownr])
  84.       Mat[i].value = -Mat[i].value;
  85.  
  86.   /* handle lower bounds by transforming to 0 */
  87.   for (i = 1; i <= Columns; i++)
  88.     if (Lowbo[Rows + i] > 0)
  89.       {
  90.     theta = Lowbo[Rows + i];
  91.     Upbo[Rows + i] -= theta;
  92.     for (j = Cend[i - 1]; j < Cend[i]; j++)
  93.       Rh[Mat[j].rownr] -= theta * Mat[j].value;
  94.       }
  95.  
  96.   for (i = 1; i <= Rows; i++)
  97.     {
  98.       num[i] = 0;
  99.       rownum[i] = 0;
  100.     }
  101.   for (i = 0; i < Nonnuls; i++)
  102.     rownum[Mat[i].rownr]++;
  103.   Rend[0] = 0;
  104.   for (i = 1; i <= Rows; i++)
  105.     Rend[i] = Rend[i - 1] + rownum[i];
  106.   for (i = 1; i <= Columns; i++)
  107.     for (j = Cend[i - 1]; j < Cend[i]; j++)
  108.       {
  109.     rownr = Mat[j].rownr;
  110.     if (rownr != 0)
  111.       {
  112.         num[rownr]++;
  113.         Colno[Rend[rownr - 1] + num[rownr]] = i;
  114.       }
  115.       }
  116.   free((char *)num);
  117.   free((char *)rownum);
  118. } /* presolve */
  119.  
  120. void  setpivcol(short *lower, 
  121.         int varin,
  122.         int numeta,
  123.         double *pcol)
  124. {
  125.   int    i, colnr;
  126.   double f;
  127.   
  128.   if (Verbose)
  129.     printf("setpivcol\n");
  130.   if ((*lower))
  131.     f = 1;
  132.   else
  133.     f = -1;
  134.   for (i = 0; i <= Rows; i++)
  135.     pcol[i] = 0;
  136.   if (varin > Rows)
  137.     {
  138.       colnr = varin - Rows;
  139.       for (i = Cend[colnr - 1]; i < Cend[colnr]; i++)
  140.     pcol[Mat[i].rownr] = Mat[i].value * f;
  141.       pcol[0] -= Extrad * f;
  142.     }
  143.   else
  144.     if ((*lower))
  145.       pcol[varin] = 1;
  146.     else
  147.       pcol[varin] = -1;
  148.   ftran(1, numeta, pcol);
  149. } /* setpivcol */
  150.  
  151.  
  152. short  colprim(int *colnr,
  153.            int *numeta,
  154.            short *minit,
  155.            double *drow)
  156. {
  157.   int    varnr, i, j;
  158.   double f, dpiv;
  159.   
  160.   if (Verbose)
  161.     printf("colprim\n");
  162.   dpiv = -EPSD;
  163.   (*colnr) = 0;
  164.   if (!(*minit))
  165.     {
  166.       for (i = 1; i <= Rows; i++)
  167.     drow[i] = 0;
  168.       drow[0] = 1;
  169.       btran((*numeta), drow);
  170.       for (i = 1; i <= Columns; i++)
  171.     {
  172.       varnr = Rows + i;
  173.       if (!Basis[varnr])
  174.         if (Upbo[varnr] > 0)
  175.           {
  176.         f = 0;
  177.         for (j = Cend[i - 1]; j < Cend[i]; j++)
  178.           f += drow[Mat[j].rownr] * Mat[j].value;
  179.         drow[varnr] = f;
  180.           }
  181.     }
  182.       for (i = 1; i <= Sum; i++)
  183.     if (abs(drow[i]) < EPSD)
  184.       drow[i] = 0;
  185.     }
  186.   for (i = 1; i <= Sum; i++)
  187.     if (!Basis[i])
  188.       if (Upbo[i] > 0)
  189.     {
  190.       if (Lower[i])
  191.         f = drow[i];
  192.       else
  193.         f = -drow[i];
  194.       if (f < dpiv)
  195.         {
  196.           dpiv = f;
  197.           (*colnr) = i;
  198.         }
  199.     }
  200.   return ((*colnr) > 0);
  201. } /* colprim */
  202.  
  203.  
  204. void  minoriteration(int colnr,
  205.              int rownr,
  206.              int *numeta)
  207. {
  208.   int    i, j, k, wk, varin, varout, elnr;
  209.   double pivot, theta;
  210.   
  211.   if (Verbose)
  212.     printf("minoriteration\n");
  213.   varin = colnr + Rows;
  214.   elnr = Endetacol[(*numeta)];
  215.   wk = elnr;
  216.   (*numeta)++;
  217.   if (Extrad != 0)
  218.     {
  219.       Eta_rownr[elnr] = 0;
  220.       Eta_value[elnr] = -Extrad;
  221.       elnr++;
  222.     }
  223.   for (j = Cend[colnr - 1] ; j < Cend[colnr]; j++)
  224.     {
  225.       k = Mat[j].rownr;
  226.       if (k == 0 && Extrad != 0)
  227.     Eta_value[Endetacol[(*numeta) - 1]] += Mat[j].value;
  228.       else
  229.     if (k != rownr)
  230.       {
  231.         Eta_rownr[elnr] = k;
  232.         Eta_value[elnr] = Mat[j].value;
  233.         elnr++;
  234.       }
  235.     else
  236.       pivot = Mat[j].value;
  237.     }
  238.   Eta_rownr[elnr] = rownr;
  239.   Eta_value[elnr] = 1 / (double) pivot;
  240.   elnr++;
  241.   theta = Rhs[rownr] / (double) pivot;
  242.   Rhs[rownr] = theta;
  243.   for (i = wk; i < elnr - 1; i++)
  244.     Rhs[Eta_rownr[i]] -= theta * Eta_value[i];
  245.   varout = Bas[rownr];
  246.   Bas[rownr] = varin;
  247.   Basis[varout] = FALSE;
  248.   Basis[varin] = TRUE;
  249.   for (i = wk; i < elnr - 1; i++)
  250.     Eta_value[i] /= - (double) pivot;
  251.   /* printf("minoriteration: new etaindex: %d\n", elnr); */
  252.   Endetacol[(*numeta)] = elnr;
  253. } /* minoriteration */
  254.  
  255.  
  256. void  rhsmincol(double *theta,
  257.         int rownr,
  258.         int varin,
  259.         int numeta)
  260. {
  261.   int    i, j, k, varout;
  262.   double f;
  263.   
  264.   if (Verbose)
  265.     printf("rhsmincol\n");
  266.   if (rownr > Rows + 1)
  267.     {
  268.       fprintf(stderr, "Error: rhsmincol called with rownr: %d, Rows: %d\n",
  269.           rownr, Rows);
  270.     }
  271.   j = Endetacol[numeta];
  272.   k = Endetacol[numeta + 1];
  273.   for (i = j; i < k; i++)
  274.     {
  275.       f = Rhs[Eta_rownr[i]] - (*theta) * Eta_value[i];
  276.       if (abs(f) < EPSB)
  277.     Rhs[Eta_rownr[i]] = 0;
  278.       else
  279.     Rhs[Eta_rownr[i]] = f;
  280.     }
  281.   Rhs[rownr] = (*theta);
  282.   varout = Bas[rownr];
  283.   Bas[rownr] = varin;
  284.   Basis[varout] = FALSE;
  285.   Basis[varin] = TRUE;
  286. } /* rhsmincol */
  287.  
  288. void  invert(int n,
  289.          int *numeta,
  290.          int *numinv)
  291. {
  292.   int    i, j, v, wk, numit, varnr, rownr, colnr, varin;
  293.   double f, theta;
  294.   double *pcol;
  295.   short  *frow;
  296.   short  *fcol;
  297.   int    *rownum, *col, *row;
  298.   int    *colnum;
  299.   short  *lbas;
  300.   double *rhs1;
  301.  
  302.   CALLOC(rownum, Rows + 1, int);
  303.   CALLOC(col, Rows + 1, int);
  304.   CALLOC(row, Rows + 1, int);
  305.   CALLOC(pcol, Rows + 2, double);
  306.   CALLOC(frow, Rows + 2, short);
  307.   CALLOC(fcol, Columns + 1, short);
  308.   CALLOC(colnum, Columns + 2, int);
  309.   CALLOC(lbas, Sum +2, short);
  310.   CALLOC(rhs1, Rows + 2, double);
  311.  
  312.   if (Verbose)
  313.     printf("a: invert:%d %d  %12f\n", (*numeta), Endetacol[(*numeta)], Rhs[0]);
  314.   
  315.   for(i = 0; i <= Sum + 1; i++)
  316.     lbas[i] = Basis[i];
  317.   for(i = 0; i <= Rows + 1; i++)
  318.     rhs1[i] = Rhs[i];
  319.   for (i = 0; i <= Rows; i++)
  320.     frow[i] = TRUE;
  321.   for (i = 0; i < n; i++)
  322.     fcol[i] = FALSE;
  323.   for (i = 0; i < Rows; i++)
  324.     rownum[i] = 0;
  325.   for (i = 1; i <= n; i++)
  326.     colnum[i] = 0;
  327.   for (i = 0; i <= Rows; i++)
  328.     if (Bas[i] > Rows)
  329.       fcol[Bas[i] - Rows - 1] = TRUE;
  330.     else
  331.       frow[Bas[i]] = FALSE;
  332.   for (i = 1; i <= Rows; i++)
  333.     if (frow[i])
  334.       for (j = Rend[i - 1] + 1; j <= Rend[i]; j++)
  335.     {
  336.       wk = Colno[j];
  337.       if (fcol[wk - 1])
  338.         {
  339.           colnum[wk]++;
  340.           rownum[i - 1]++;
  341.         }
  342.     }
  343.   for (i = 1; i <= Rows; i++)
  344.     Bas[i] = i;
  345.   for (i = 1; i <= Rows; i++)
  346.     Basis[i] = TRUE;
  347.   for (i = 1; i <= n; i++)
  348.     Basis[i + Rows] = FALSE;
  349.   for (i = 0; i <= Rows; i++)
  350.     Rhs[i] = Rh[i];
  351.   for (i = 1; i <= n; i++)
  352.     {
  353.       varnr = Rows + i;
  354.       if (!Lower[varnr])
  355.     {
  356.       theta = Upbo[varnr];
  357.       for (j = Cend[i - 1]; j < Cend[i]; j++)
  358.         Rhs[Mat[j].rownr] -= theta * Mat[j].value;
  359.     }
  360.     }
  361.   for (i = 1; i <= Rows; i++)
  362.     if (!Lower[i])
  363.       Rhs[i] -= Upbo[i];
  364.   (*numeta) = 0;
  365.   v = 0;
  366.   rownr = 0;
  367.   (*numinv) = 0;
  368.   numit = 0;
  369.   while (v < Rows)
  370.     {
  371.       rownr++;
  372.       if (rownr > Rows)
  373.     rownr = 1;
  374.       v++;
  375.       if (rownum[rownr - 1] == 1)
  376.     if (frow[rownr])
  377.       {
  378.         v = 0;
  379.         j = Rend[rownr - 1] + 1;
  380.         while (!(fcol[Colno[j] - 1]))
  381.           j++;
  382.         colnr = Colno[j];
  383.         fcol[colnr - 1] = FALSE;
  384.         colnum[colnr] = 0;
  385.         for (j = Cend[colnr - 1]; j < Cend[colnr]; j++)
  386.           if (frow[Mat[j].rownr])
  387.         rownum[Mat[j].rownr - 1]--;
  388.         frow[rownr] = FALSE;
  389.         minoriteration(colnr, rownr, numeta);
  390.       }
  391.     }
  392.   v = 0;
  393.   colnr = 0;
  394.   while (v < n)
  395.     {
  396.       colnr++;
  397.       if (colnr > n)
  398.     colnr = 1;
  399.       v++;
  400.       if (colnum[colnr] == 1)
  401.     if (fcol[colnr - 1])
  402.       {
  403.         v = 0;
  404.         j = Cend[colnr - 1] + 1;
  405.         while (!(frow[Mat[j - 1].rownr]))
  406.           j++;
  407.         rownr = Mat[j - 1].rownr;
  408.         frow[rownr] = FALSE;
  409.         rownum[rownr - 1] = 0;
  410.         for (j = Rend[rownr - 1] + 1; j <= Rend[rownr]; j++)
  411.           if (fcol[Colno[j] - 1])
  412.         colnum[Colno[j]]--;
  413.         fcol[colnr - 1] = FALSE;
  414.         numit++;
  415.         col[numit - 1] = colnr;
  416.         row[numit - 1] = rownr;
  417.       }
  418.     }
  419.   for (j = 1; j <= n; j++)
  420.     if (fcol[j - 1])
  421.       {
  422.     fcol[j - 1] = FALSE;
  423.     setpivcol(&Lower[Rows + j], j + Rows, (*numeta), pcol);
  424.     rownr = 1;
  425.     while (!(frow[rownr] && pcol[rownr]))
  426.       rownr++; /* this sometimes sets rownr to Rows + 2 and makes */
  427.            /* rhsmincol crash. MB */
  428.     frow[rownr] = FALSE;
  429.     condensecol(rownr, (*numeta), pcol);
  430.     theta = Rhs[rownr] / (double) pcol[rownr];
  431.     rhsmincol(&theta, rownr, Rows + j, (*numeta));
  432.     addetacol(numeta);
  433.       }
  434.   for (i = numit - 1; i >= 0; i--)
  435.     {
  436.       colnr = col[i];
  437.       rownr = row[i];
  438.       varin = colnr + Rows;
  439.       for (j = 0; j <= Rows; j++)
  440.     pcol[j] = 0;
  441.       for (j = Cend[colnr - 1]; j < Cend[colnr]; j++)
  442.     pcol[Mat[j].rownr] = Mat[j].value;
  443.       pcol[0] -= Extrad;
  444.       condensecol(rownr, (*numeta), pcol);
  445.       theta = Rhs[rownr] / (double) pcol[rownr];
  446.       rhsmincol(&theta, rownr, varin, (*numeta));
  447.       addetacol(numeta);
  448.     }
  449.   for (i = 1; i <= Rows; i++)
  450.     if (abs(Rhs[i]) < EPSB)
  451.       Rhs[i] = 0;
  452.   f = 0;
  453.   for (i = 1; i <= Rows; i++)
  454.     if (Rhs[i] < 0)
  455.       f += Rhs[i];
  456.     else
  457.       if (Rhs[i] > Upbo[Bas[i]])
  458.         f = f + Upbo[Bas[i]] - Rhs[i];
  459.   if (Verbose)
  460.     printf("b: invert:%d %d  %12f  %12f\n", (*numeta), Endetacol[(*numeta)],
  461.        Rhs[0], f);
  462.  
  463.   free((char *)rownum);
  464.   free((char *)col);
  465.   free((char *)row);
  466.   free((char *)pcol);
  467.   free((char *)frow);
  468.   free((char *)fcol);
  469.   free((char *)colnum);
  470.   free((char *)lbas);
  471.   free((char *)rhs1);
  472. } /* invert */
  473.  
  474.  
  475. short  rowprim(int *rownr,
  476.            double *theta,
  477.            double *pcol)
  478. {
  479.   int    i;
  480.   double f, quot;
  481.   
  482.   if (Verbose)
  483.     printf("rowprim\n");
  484.   (*rownr) = 0;
  485.   (*theta) = INFINITE;
  486.   for (i = 1; i <= Rows; i++)
  487.     {
  488.       f = pcol[i];
  489.       if (f != 0)
  490.     {
  491.       quot = (*theta) * 10;
  492.       if (f > 0)
  493.         quot = Rhs[i] / (double) f;
  494.       else
  495.         if (Upbo[Bas[i]] < INFINITE)
  496.           quot = (Rhs[i] - Upbo[Bas[i]]) / (double) f;
  497.       if (quot < (*theta))
  498.         {
  499.           (*theta) = quot;
  500.           (*rownr) = i;
  501.         }
  502.     }
  503.     }
  504.   return ((*rownr) > 0);
  505. } /* rowprim */
  506.  
  507.  
  508.  
  509. void  iteration(int *numc,
  510.         int *rownr,
  511.         int *varin,
  512.         double *theta,
  513.         double *up,
  514.         short *minit,
  515.         short *low,
  516.         short *primair,
  517.         short *smotes,
  518.         int *iter,
  519.         int *numinv)
  520. {
  521.   int    i, k, varout;
  522.   double f, pivot;
  523.   
  524.   if (Verbose)
  525.     printf("iteration\n");
  526.   (*iter)++;
  527.   (*minit) = (*theta) > (*up);
  528.   if ((*minit))
  529.     {
  530.       (*theta) = (*up);
  531.       (*low) = !(*low);
  532.     }
  533.   k = Endetacol[(*numc) + 1];
  534.   pivot = Eta_value[k - 1];
  535.   for (i = Endetacol[(*numc)]; i < k; i++)
  536.     {
  537.       f = Rhs[Eta_rownr[i]] - (*theta) * Eta_value[i];
  538.       if (abs(f) < EPSB)
  539.     Rhs[Eta_rownr[i]] = 0;
  540.       else
  541.     Rhs[Eta_rownr[i]] = f;
  542.     }
  543.   if (!(*minit))
  544.     {
  545.       Rhs[(*rownr)] = (*theta);
  546.       varout = Bas[(*rownr)];
  547.       Bas[(*rownr)] = (*varin);
  548.       Basis[varout] = FALSE;
  549.       Basis[(*varin)] = TRUE;
  550.       if ((*primair) && pivot < 0)
  551.     Lower[varout] = FALSE;
  552.       if (!(*low) && (*up) != INFINITE)
  553.     {
  554.       (*low) = TRUE;
  555.       Rhs[(*rownr)] = (*up) - Rhs[(*rownr)];
  556.       for (i = Endetacol[(*numc)]; i < k; i++)
  557.         Eta_value[i] = -Eta_value[i];
  558.     }
  559.       addetacol(numc);
  560.       (*numinv)++;
  561.     }
  562.   if ((*smotes) && Verbose)
  563.     {
  564.       printf("Iteration %d: ", (*iter));
  565.       if ((*minit))
  566.     {
  567.       printf("%4d", (*varin));
  568.       if (Lower[(*varin)])
  569.         printf("  ");
  570.       else
  571.         printf(" u");
  572.       printf("%22c", ' ');
  573.     }
  574.       else
  575.     {
  576.       printf("%4d", (*varin));
  577.       if (Lower[(*varin)])
  578.         printf("  ");
  579.       else
  580.         printf(" u");
  581.       printf("%4d", varout);
  582.       if (Lower[varout])
  583.         printf("  ");
  584.       else
  585.         printf(" u");
  586.       printf("%12.5f", pivot);
  587.     }
  588.       printf("    %12f", Pcol[0]);
  589.       if (!(*primair))
  590.     {
  591.       f = 0;
  592.       for (i = 1; i <= Rows; i++)
  593.         if (Rhs[i] < 0)
  594.           f += Rhs[i];
  595.         else
  596.           if (Rhs[i] > Upbo[Bas[i]])
  597.         f = f + Upbo[Bas[i]] - Rhs[i];
  598.       printf("  %12f", f);
  599.     }
  600.       else
  601.     printf("    %12f", Rhs[0]);
  602.       printf("\n");
  603.     }
  604. } /* iteration */
  605.  
  606.  
  607. int  solvelp(void)
  608. {
  609.   int    i, j, iter, varnr;
  610.   double f, theta;
  611.   short  primair, doiter, doinvert, smotes;
  612.   double *drow, *prow;
  613.   int    numinv, invnum;
  614.   int    numeta;
  615.   short  artif, minit;
  616.   int    colnr, rownr;
  617.  
  618.   CALLOC(drow, Sum + 2, double);
  619.   CALLOC(prow, Sum + 2, double);
  620.  
  621.   numinv = 0;
  622.   invnum = 50; /* number of iterations between inversions */
  623.   numeta = 0;
  624.   iter = 0;
  625.   minit = FALSE;
  626.   primair = TRUE;
  627.   i = 0;
  628.   smotes = TRUE;
  629.  
  630.   while (i != Rows && primair)
  631.     {
  632.       i++;
  633.       primair = Rhs[i] >= 0 && Rhs[i] <= Upbo[Bas[i]];
  634.     }
  635.   
  636.   if (!primair)
  637.     {
  638.       drow[0] = 1;
  639.       for (i = 1; i <= Rows; i++)
  640.     drow[i] = 0;
  641.       Extrad = 0;
  642.       for (i = 1; i <= Columns; i++)
  643.     {
  644.       varnr = Rows + i;
  645.       drow[varnr] = 0;
  646.       for (j = Cend[i - 1]; j < Cend[i]; j++)
  647.         if (drow[Mat[j].rownr] != 0)
  648.           drow[varnr] += drow[Mat[j].rownr] * Mat[j].value;
  649.       if (drow[varnr] < Extrad)
  650.         Extrad = drow[varnr];
  651.     }
  652.       artif = Extrad < -EPSD;
  653.     }
  654.   else
  655.     {
  656.       artif = FALSE;
  657.       Extrad = 0;
  658.     }
  659.   if (Verbose)
  660.     printf("artificial:%12.5f\n", Extrad);
  661.   minit = FALSE;
  662.   do {
  663.     doiter = FALSE;
  664.     doinvert = FALSE;
  665.     if (primair)
  666.       {
  667.     rownr = 0;
  668.     if (colprim(&colnr, &numeta, &minit, drow))
  669.       {
  670.         setpivcol(&Lower[colnr], colnr, numeta, Pcol);
  671.         if (rowprim(&rownr, &theta, Pcol))
  672.           {
  673.         doiter = TRUE;
  674.         condensecol(rownr, numeta, Pcol);
  675.           }
  676.       }
  677.       }
  678.     else
  679.       {
  680.     if (!minit)
  681.       rowdual(&rownr);
  682.     if (!doinvert)
  683.       if (rownr > 0)
  684.         {
  685.           if (coldual(&numeta, &rownr, &colnr, &minit, prow, drow))
  686.         {
  687.           doiter = TRUE;
  688.           setpivcol(&Lower[colnr], colnr, numeta, Pcol);
  689.           condensecol(rownr, numeta, Pcol);
  690.           f = Rhs[rownr] - Upbo[Bas[rownr]];
  691.           if (f > 0)
  692.             {
  693.               theta = f / (double) Pcol[rownr];
  694.               if (theta <= Upbo[colnr])
  695.             Lower[Bas[rownr]] = !Lower[Bas[rownr]];
  696.             }
  697.           else /* getting div by zero here ... MB */
  698.             theta = Rhs[rownr] / (double) Pcol[rownr];
  699.         }
  700.         }
  701.       else
  702.         if (artif)
  703.           {
  704.         primair = TRUE;
  705.         artif = FALSE;
  706.         printf("primal\n");
  707.         doinvert = TRUE;
  708.         Extrad = 0;
  709.           }
  710.         else
  711.           colnr = 0;
  712.       }
  713.     if (doiter)
  714.       iteration(&numeta, &rownr, &colnr, &theta, &Upbo[colnr],
  715.         &minit, &Lower[colnr], &primair, &smotes, &iter, &numinv);
  716.     if (numinv >= invnum)
  717.       doinvert = TRUE;
  718.     if (doinvert)
  719.       invert(Columns, &numeta, &numinv);
  720.   } while (rownr && colnr || doinvert);
  721.  
  722.   free((char *)drow);
  723.   free((char *)prow);
  724.  
  725.   return(rownr || colnr); /* return nonzero if anything was wrong */
  726. } /* solvelp */
  727.  
  728.  
  729. void construct_solution(double *sol)
  730. {
  731.   int    i, j, basi;
  732.   double f;
  733.  
  734.   /* zero all results of rows */
  735.   bzero((char *)sol, (Rows + 1) * sizeof(double));
  736.  
  737.   for (i = Rows + 1; i <= Sum; i++)
  738.     sol[i] = Lowbo[i];
  739.   for (i = 1; i <= Rows; i++)
  740.     {
  741.       basi = Bas[i];
  742.       if (basi > Rows)
  743.     sol[basi] += Rhs[i];
  744.     }
  745.   for (i = Rows + 1; i <= Sum; i++)
  746.     if (!Basis[i] && !Lower[i])
  747.       sol[i] += Upbo[i];
  748.   for (j = 1; j <= Columns; j++)
  749.     {
  750.       f = sol[Rows + j];
  751.       if (f != 0)
  752.     for (i = Cend[j - 1]; i < Cend[j]; i++)
  753.       sol[Mat[i].rownr] += f * Mat[i].value;
  754.     }
  755.   
  756.   for (i = 0; i <= Rows; i++)
  757.     {
  758.       if (abs(sol[i]) < EPSB)
  759.     sol[i] = 0;
  760.       else
  761.     if (Chsign[i])
  762.       sol[i] = -sol[i];
  763.     }
  764. } /* construct_solution */
  765.  
  766.  
  767. int solve(double *upbo,
  768.       double *lowbo)
  769. {
  770.   int    i, failure, notint;
  771.   intrec *ptr;
  772.  
  773.   Level++;
  774.   debug_print("starting solve");
  775.  
  776.   /* make fresh copies of Mat, Rh, Upbo, Lowbo, as solving changes them */
  777.   bcopy((char*)Orig_mat,  (char *)Mat,   (Nonnuls + 1) * sizeof(matrec));
  778.   bcopy((char*)Orig_rh,   (char *)Rh,    (Rows + 2)    * sizeof(double));
  779.   bcopy((char*)upbo,      (char *)Upbo,  (Sum + 2)     * sizeof(double));
  780.   bcopy((char*)lowbo,     (char *)Lowbo, (Sum + 2)     * sizeof(double));
  781.  
  782.   Rend[0] = 0;
  783.   presolve();
  784.  
  785.   for (i = 0; i <= Rows; i++)
  786.     Rhs[i] = Rh[i];
  787.  
  788.   Endetacol[0] = 0;
  789.   failure = solvelp();
  790.  
  791.   if(failure)
  792.     debug_print("this problem has no solution");
  793.  
  794.   if(!failure) /* there is a solution */
  795.     {
  796.       construct_solution(Solution);
  797.  
  798.       debug_print("a solution was found");
  799.       debug_print_solution(Solution);
  800.  
  801.       /* if this solution is worse than the best sofar, this branch must die */
  802.       if(Solution[0] <= Best_solution[0])
  803.     {
  804.       debug_print("but it was worse than the best sofar, discarded");
  805.       Level--;
  806.       return(TRUE);
  807.     }
  808.  
  809.       /* check if solution contains enough ints */
  810.       notint = 0;
  811.       for(ptr = First_int; !notint && ptr; ptr = ptr->next)
  812.     if(!is_int(Solution[ptr->varnr]))
  813.       notint = ptr->varnr;
  814.  
  815.       if(notint) /* there is at least one value not yet int */
  816.     {
  817.       /* set up two new problems */
  818.       double *new_upbo, *new_lowbo;
  819.       double new_bound;
  820.       int    res1, res2;
  821.  
  822.       /* allocate room for them */
  823.       MALLOC(new_upbo,  Sum + 2, double);
  824.       MALLOC(new_lowbo, Sum + 2, double);
  825.  
  826.       bcopy((char *)upbo,  (char *)new_upbo,  (Sum + 2) * sizeof(double));
  827.       bcopy((char *)lowbo, (char *)new_lowbo, (Sum + 2) * sizeof(double));
  828.  
  829.       debug_print("not enough ints. Selecting var %s, val: %10.3g",
  830.               Names[notint], Solution[notint]);
  831.  
  832.       new_bound = floor(Solution[notint]);
  833.       new_upbo[notint] = new_bound;
  834.  
  835.       debug_print("starting first subproblem with bounds:");
  836.       debug_print_bounds(new_upbo, lowbo);
  837.  
  838.       res1 = solve(new_upbo, lowbo);
  839.  
  840.       new_bound += 1;
  841.       new_lowbo[notint] = new_bound;
  842.  
  843.       debug_print("starting second subproblem with bounds:");
  844.       debug_print_bounds(upbo, new_lowbo);
  845.  
  846.       res2 = solve(upbo, new_lowbo);
  847.  
  848.       if(res1 && res2) /* both failed */
  849.         failure = TRUE;
  850.       else
  851.         failure = FALSE;
  852.  
  853.       FREE(new_upbo);
  854.       FREE(new_lowbo);
  855.     }
  856.       else /* all required values are int */
  857.     {
  858.       debug_print("--> Valid solution found");
  859.       if(Solution[0] > Best_solution[0]) /* Current solution better */
  860.         {
  861.           debug_print("the best sofar. Prev: %10.3g, New: %10.3g",
  862.               Best_solution[0], Solution[0]);
  863.           bcopy((char *)Solution, (char *)Best_solution,
  864.             (Sum + 1) * sizeof(double));
  865.           if(Show_results)
  866.         {
  867.           fprintf(stderr, "Intermediate solution:\n");
  868.           print_solution(stderr, Best_solution);
  869.         }
  870.         }
  871.     }
  872.     }
  873.   Level--;
  874.   return(failure);
  875. } /* solve */
  876.